Accurate sampling using Langevin dynamics 
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We show how to derive a simple integrator for the Langevin equation and illustrate how it is 
possible to check the accuracy of the obtained distribution on the fly, using the concept of effective 
energy introduced in a recent paper [J. Chem. Phys. 126, 014101 (2007)]. Our integrator leads 
to correct sampling also in the difficult high-friction limit. We also show how these ideas can be 
applied in practical simulations, using a Lennard- Jones crystal as a paradigmatic case. 



I. INTRODUCTION 



Langevin dynamics was first introduced in molecular 
simulations to calculate the properties of mesoscopic sys- 
tems [l[ . Here a dissipative force and a noise were added 
to the Hamilton equations to model a bath of lighter 
particles. The formal justification for this model can be 
obtained using the projection operator techniques @, Q. 
However, it was soon realized that Langevin dynamics 
can also be used as a thermostat [H, adding the dissipa- 
tive forces and the noise to the Hamiltonian dynamics to 
allow a molecular dynamics simulation to explore an en- 
semble at a fixed temperature. Furthermore, it has been 
used to sample arbitrary distribution, for instance in the 
case of numerical quantum-chromodynamics 5] . 

Several algorithms have been proposed for the numeri- 
cal integration of the Langevin equation, see among oth- 
ers Refs. SSSStMEtilMQIilMEtili. Most 

of them were derived with the aim of producing accurate 
trajectories, i.e. dynamical properties, up to a given or- 
der. Because of that, they usually break down when a 
high friction is applied, essentially when the velocities 
are varying too fast with respect to the chosen time step. 
Moreover, their design is not focused on the correctness 
of the ensemble generated. A notable exception is given 
by the schemes derived in Ref. (l7| , where the free param- 
eters of the algorithm are chosen so as to minimize the 
sampling errors. However, none of the algorithms so far 
proposed offer any way of checking the accuracy of the 
sampling during a numerical simulation. This is at vari- 
ance with the numerical integration of Hamilton's equa- 
tions, where the conservation of the total energy has been 
traditionally used to this end [T^, Gj[ . The standard ap- 
proach in molecular dynamics is thus to choose the time 
step by monitoring the energy conservation in a few mi- 
crocanonical runs, then to adopt the same time step for 
the Langevin dynamics. To the best of our knowledge, 
only in a recent paper [i~5| Scemama et al. have shown 
how to correct exactly the discretization errors in the 
Langevin dynamics in the context of variational Monte 
Carlo, using a Metropolis procedure. However, the poor 



scaling of these accept-reject algorithms with respect to 
the number of degrees of freedom prevents thei r ap plica- 
tion to global moves in very large systems [l^, l2d| |. 

In a recent paper [2l[ we introduced a constant- 
temperature molecular-dynamics method. In that con- 
text, we discussed the notion of effective energy, as a 
measure of sampling accuracy. In Ref. (2lTj only one vari- 
able, the total kinetic energy, was subject to stochastic 
fluctuations and the response of the thermostat could be 
modeled so as to have a minimal effect on the dynamics. 
Here we apply some of the ideas developed in Ref. [2l| to 
Langevin dynamics, where all degrees of freedom can be 
separately controlled and the time scale over which the 
thermostat reacts is defined by the friction coefficient. 
When used as thermostat, Langevin dynamics can be 
more efficient in difficult cases, but it is more disruptive 
of the dynamics. In an extension of Ref. [2l[, we integrate 
Langevin using a simple algorithm derived from a Trot- 
ter decomposition. The effective-energy drift allows the 
sampling error to be controlled during a simulation, and 
can be used in a rigorous way to perform reweighting or 
accept-reject algorithms, in a scheme that turns out to be 
similar to that discussed by Scemama et al. [IB] . The ad- 
vantage of our formulation is that, for large systems, the 
effective energy can be simply checked against long-term 
drifts, in the same way as the total energy has tradition- 
ally been used to check the accuracy of microcanonical 
molecular dynamics. We also show the properties of the 
effective energy in model harmonic oscillators and in a 
realistic Lennard- Jones crystal. 



II. THEORY 
A. Langevin dynamics 

We consider a particle with mass m subject to a poten- 
tial energy U(q). The generalization to multiple degrees 
of freedom is straightforward. The probability density 
for the canonical ensemble at an inverse temperature j3 
is 
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P(p,q) dpdpot e -P^e~P u{q) dpdp. 



(1) 
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The canonical ensemble can be sampled through the 
Langevin dynamics 



*(*)=/(«(*)) dt-jp(t) dt 

dq(t) = dt, 
m 



2mj 



dW(t) 



(2a) 
(2b) 



where f(q) = is the deterministic force, 7 is the fric- 
tion coefficient, and dW(t) is a Wiener noise in the Itoh 
convention [Hj], normalized as (dW(t)dW(t')) = S(t-t'). 
A description equivalent to the stochastic Eq. (|2|) can 
be formulated in terms of the probability density, which 
evolves according to the Fokker-Planck equation [22, HH 



which are defined as 
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dp 
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(8a) 
(8b) 
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Several choices are now available for the Trotter splitting. 
We notice that the operators e~^ h ~< and e ~ AtL pi leave 
the stationary distribution in Eq. (TTJ) unchanged: 



- — L_ 



P = P ; e 



-AtL„ 



P = P. 



(9) 



dP(p,g;t) 
dt 



LP(p,q;t) 



where 



(3) 



t ,/ % d p d ( d ra<9, 
L = f{q) Yp + V l dq-- 1 \d-p P+ ^W ] - () 



This is due to the fact that the canonical distribution 
is stationary not only with respect to L but also with 
respect to L pq = L p + L q , which corresponds to Hamilton 
propagation, and with respect to L 7 , which introduces 
the combined effect of friction and noise. Thus, even 
if the commutator [L pqi Lj] ^ 0, the following splitting 
does not introduce sampling errors, 



The formal solution of Eq. © at a finite time step At is 



-AtL 



At T A + T At J 



(10) 



P(p,q;t + At) 



-AtL 



P(p,q;t), 



(5) 



which however cannot be evaluated explicitly. Notice 
that for Hamiltonian dynamics, 7 = 0, the operator 
L is anti-Hermitian and the propagator e 



-AtL 



is uni- 



tary. These properties hold only for a deterministic area- 
preserving dynamics. They do not hold in a Langevin 
process. 



B. A simple integrator 

As was first recognized by Tuckerman et al. (24j and, 
independently, by Sexton and Weingarten (25[, the Trot- 
ter formula [261 ] allows an approximated propagator to be 
constructed as 



M 



-AtL 



n ^ n 



At T 

, o~-kfc 



(6) 



-M 



fc=l 



where M is the number of stages in the integrator and 
^2 j Lj = L. Since in general the Lj's do not commute 
among themselves, the order in which the stages are ap- 
plied is relevant, and the splitting in Eq. ([6]) introduces 
some error into the propagation. The key point here is 
that the stages e~~ Lj are chosen so that they can be 
integrated analytically, and the Trotter splitting is the 
only source of errors. 

It is natural to write L as a sum of three parts: 



(7) 



since it can be interpreted as a sequence of moves each of 
which has the correct limiting distribution. The e~~ L '< 
move provides ergodicity in the momenta subspace only, 
while the e~ AtLpq move mixes the momenta and positions 
subspaces. An integrator designed to apply the propaga- 
tor in Eq. IjlOp would provide an approximate trajectory 
and an exact sampling, independently of At and 7. The 
propagator e~^~ L ~< can be integrated analytically. Un- 
fortunately, the propagator e - AtL Pi cannot be integrated 
exactly and has to be split further. We opt here for the 
simplest choice, which is the same used to obtain the 
velocity Verlet algorithm: 



-AtL 



-A±T 



-AtL„ 



(11) 



In specific cases, different decompositions of L pq could 
be adopted. For example, if the forces can be separated 
into contributions varying on different time scales, a 
multiple-time-step decomposition is expected to be more 
efficient 0]. 

Other possible choices for the Trotter splitting which 
are substantially equivalent to Eq. (fTTj) can be obtained, 



based on the three operators L q , L p 



and L. 



It 



worthwhile to notice that in principle there is no need 



tO Split u p aim jj 7 , 



and , since L p -^ — L p 



L 7 can be 



also evolved analytically. Ricci and Ciccotti [14( de- 
rived two integrators using splittings that, in our nota- 



tion, would read e 

-AtL 



-AtL 



-AtL „ 



At T 
-1 T-J-'O 



"e ""-^se 2 L Pt and 
e ps e A r L i e - AtL Pi e ^ A r L q . These decompositions 
involve a single splitting and thus appear more accurate 
than Eq. (fTTj) . However, when 7 At is negligible, they do 
not offer any advantage, and when 7At is not negligible, 



3 



they do not sample the proper ensemble. This can be 
easily verified taking the limit 7 At — + oo. On the other 
hand, in our scheme the only ensemble violations arise 
from the fact that for a finite At the evolution of L pq is 
approximated. These violations are independent of the 
choice of the friction. Even the infinite friction limit can 
be taken safely, as shown in Appendix[XJ Thus, when the 
sampling quality is an issue, our scheme offers significant 
advantages. 

The splitting in Eq. (fTTj) leads to an explicit integration 
scheme. In the derivation we use the analytical propa- 
gation formula for L 1 which can be found in Ref. [231 ] . 
After some manipulation, the integrator is written as: 



P {t+) =c lP (t)+c 2 R(t) 

q(t + At) = q(t) + ^-At 



p{t~ + At) =p{t+) + 



772 777 2 

/(«(*)) + /(«(* + At)) 



p(t + At) = cip{r + At) + c 2 R'{t + At) 



(12a) 
(12b) 

At (12c) 
(12d) 



where R and R' are two independent Gaussian numbers 
and the coefficients c\ and c 2 are 



Cl 



C2 







(13a) 



(13b) 



Equation (|13b|) fixes the weight of the rescaling factor c\ 
and of the amplitude of the Gaussian number c 2 in such 
a way that C\p + c 2 i? will be distributed in the same way 
as p. Thus, Eq. (|13b[) alone guarantees the correctness 
of the sampling. On the other hand, Eq. (|13a|) gives the 
relation between the friction 7 and the rescaling factor 

Cl. 

In Equation (|12p. the combination of the two inner 
stages is a velocity Verlet step, and corresponds to the 
approximate propagation of e~ AtLpq . The first and last 
stages represent the action of the thermostat, i.e. the 
exact propagation of e~~ tL ~< . We denote as p(t + ) and 
p(t~) the momenta immediately after and immediately 
before the action of the thermostat. We also observe that 
the first and last stages can be merged as p(t + + At) = 
c?p(t~ + At) + c 2 a/c? + lR(t + At) so that one Gaus- 
sian random number per degree of freedom is required 
at each step. This allows the simulation to speed up 
when the calculation of the deterministic forces is partic- 
ularly cheap and the generation of the Gaussian random 
numbers becomes computationally relevant. If one is in- 
terested in the values of the momenta at time t, i.e. syn- 
chronized with the positions, they can be reconstructed 
afterwards. 



C. Control of sampling errors 

We now use the concept of effective energy H intro- 
duced in Ref. 21] to control the accuracy of the sam- 



pling. For clarity we repeat here some of the notions 
already presented there. 

Our goal is to generate a sequence of points Xi = (j>i, qi) 
in the phase-space, so that a time average can be used 
in place of the ensemble average 0. Usually, in molec- 
ular dynamics simulations this sampling is approximate, 
due to the finite-time-step errors. On the other hand, 
in a Monte Carlo simulation the moves are accepted or 
refused in such a way that the exact distribution is en- 
forced. Here, we interpret a stochastic molecular dy- 
namics as a highly efficient Monte Carlo where all the 
moves are accepted. We define M(a;,+i «— X4) dxi+\ 
the distribution probability of the point OJj+i to be cho- 
sen as the next point, given that the present point is x^. 
We also define the conjugate point x* — (~p,q), which 
is obtained by inverting the momentum, and satisfies 
P(x) = P(x*). If EquationJ© was integrated exactly, 
then the detailed balance [22| would be satisfied, i.e., 
M(x i+1 «- Xi)P(Xi) = M(x* «- x* +1 )P(x* +1 ). However, 
this is not true when a finite time step is used. Thus, we 
introduce a weight u>j associated to the point Xi, which 
evolves as 



w i+ i _ M{x* «- x* l+1 )P{x l+l ) 



M(x 



Xi)P(Xi) 



(14) 



The same information can be expressed in terms of an 
effective energy, defined as Hi = — i log Wi , which evolves 
according to 



H i+1 -Hi = --log ( -^jj- )+H(x i+1 )-H(xi 



13 °\M{x i+ i^Xi) 

(15) 

We now proceed into an explicit derivation of the terms 
needed. 

In standard hybrid Monte Carlo, the trial moves 
are generated using an area-preserving scheme, so that 
M(x* <— x* +l ) = M(x l+1 <— Xi). Thus, the effective 
energy H reduces to the Hamiltonian H. However, the 
Langevin equation is explicitly non-area-preserving, and 
an additional contribution due to phase-space compres- 
sion has to be evaluated. We now calculate it explicitly 
for the integrator in Eq. (fT2")) . In Ref. [2l[ we used the 
fact that the thermostat moves are designed so as to sat- 
isfy detailed balance. We present here a more general 
way of evaluating this contribution that can be straight- 
forwardly applied to other integrators. 

We recall that the random numbers R and R' are 
drawn from a Gaussian distribution, i.e. 



P(R, R') dRdR' 



1 

27 



— dRdR'. 



(16) 



We notice that given the starting point x^ — (pi,qi) and 
the ending point = {pi+i, qi+i) the value of R and 
R' can be determined solving Eqs. (|12p with respect to 
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R and R': 



R 



R' 



(ft+i - 
-(ft+i 



ft) 



c 2 At 
C\m 
c 2 At 



/(gi)At ci 

— ^ Pi 

_ ci/(g i+ i)Af 
2c 2 



C2 



"Pi 



(17a) 
(17b) 



where we have identified the sequence index i with the 
time t and the sequence index i + 1 with the time t + Ai. 
Now, changing the variables from (R,R') to (gj_|_i,pj_)_x) 
one obtains the following expression for the transition 
probability: 



AI 



\{Qi+i,Pi+i) <- (ft, Pi)) 



m / 1 / rn /fa) At ^ 2 

^r^il^-^A* — 2 — cm 



2itc\At 



^2 ((ft+1 - ft) 



"aT 



f(q i+1 ) Cl At 

o Pt+1 



(18) 



At this stage we know the probability for the forward 
move M((qi + i,pi + i) <— (ft, Pi)). With a similar proce- 
dure we can find the probability for the backward move, 
M{{qi,—pi) <— , — pi+i)), and, with some further 

manipulation, the contribution of the phase-space com- 
pression to the effective energy: 



1 ( M(x* <- 
13 ° S \M(x i+1 



X i+l) 



+ (ft+1 - ft) 



Xi) 
fi+1 + fi 



1? 

Pi+1 



2m 



Pi 
2m 



At 2 



(/(<?■ 



i+i; 



/(ft) 2 )- 



(19) 



For this derivation it is crucial that the change of vari- 
ables be well defined. Since we have two noise terms 
(R, R') and two variables (qi+i,pi+\), we have to require 
the Jacobian of the transformation to be different from 
zero. For integrators using only one noise term, it is not 
obvious that, given the forward trajectory, the backward 
trajectory is possible. If the backward trajectory is pos- 
sible, then the effective-energy drift depends on the ra- 
tio between the forward and backward probabilities and 
gives a quantitative measure of the violation of detailed 
balance. If the backward trajectory is not possible, then 
the integrator cannot satisfy detailed balance. As an ex- 
ample, the second integrator introduced by Ricci and Ci- 
ccotti cannot satisfy detailed balance, as was already 
pointed out by Scemama et al. [laL On the other hand, 
the modification described in Ref. [15| can satisfy detailed 
balance. It is interesting that in Ref. [l5| the authors are 
using the usual formulation of detailed balance, which 
leads to the need for an explicit inversion of the sign of 
the velocities. We use a more general formulation of de- 
tailed balance [22| in which velocities are considered as 
odd variables and their inversion after an accepted step 
is not required. One could also object that the condition 



of detailed balance is not strictly necessary [27j. How- 
ever, it appears to us that detailed balance is the only 
way to enforce or check a distribution in a local manner, 
i.e., using only information about the present point Xi, 
the next point Xi+i, and their conjugated points x* and 

X i+V 

Equations (| 1 5[) and (flT))) can be combined, giving a final 
expression for the effective-energy increment as 



AH 



Aq{^ 



i+l, 



At 2 

+ AU+— A(/ 2 ). (20) 
8m 



From Equation (|20[) it is easy to see that when At is small 
enough the effective energy is approximately constant, 
since the first and second terms tend to compensate each 
other and the third term vanishes on the order of At 2 . 
We also notice that the third term in Eq. (f2"0")) is an exact 
differential. Thus it contributes to the fluctuations of the 
effective energy but not to its drift. 

We notice that the increment of the effective energy 
in Eq. Ip20)) is exactly equal to the difference of the total 
energy before and after the velocity Verlet step, as in the 
case of the scaling procedure described in Ref. [2l|. In 
fact, also here we can think of our dynamics as composed 
of a combination of two steps: one, described by the op- 
erator e~^ L i , which exactly satisfied detailed balance; 
the other, which is the velocity Verlet step, is symplec- 
tic but does not exactly conserve the energy. Only the 
violations arising from the latter are accumulated into 
the effective energy. Thus, in practice, H can be cal- 
culated simply by summing the increments of the total 
energy due to the Verlet, discarding the increments due 
to the thermostat. Alternatively, it can be obtained by 
subtracting from the total energy the sum of all its incre- 
ments due to the thermostat. The same procedure can be 
applied directly also in the case of the Peters thermostat 
[28j |. based on dissipative particle dynamics [29l |. where 
e -^r L t i s substituted by a rescaling of the relative ve- 
locity of neighboring particles, the only condition being 
the fact that the rescalings are performed in a way that 
analytically preserves the target ensemble. 

The effective energy can be calculated on the fly and, 
aside from numerical truncation errors, it gives a quan- 
titative way to assess the accuracy of the calculation. In 
the spirit of Ref. [3fl | , one can obtain exact ensemble av- 
erages with (A) — J2i w iA(xi)/J2i w i- The variation of 
H on segments of trajectory can also be used in a hy- 
brid Monte Carlo scheme [20|], where the acceptance is 

calculated as min 

f l e -0A(H)\ j n thig latter 

case, our 

scheme becomes similar to that presented by Scemama 
et al. [l5l |. In a molecular dynamics context, the effective 
energy H is simply monitored during the simulation. It 
may fluctuate, but it should not exhibit a large system- 
atic drift. 
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III. EXAMPLES 



A. Harmonic oscillator 



It is instructive to study the properties of the integra- 
tor in Eq. (|12[) when it is applied to a harmonic oscillator. 
We consider an energy profile 



U{q) = imwV- 



(21) 



We are interested in the time evolution of the effective 
energy H. It can be easily shown that for a quadratic 
potential the first two terms in Eq. (|20|) cancel exactly, 
and only the third term survives. Thus, the integral over 
the trajectory is not necessary and the effective energy 
H is a state function 



At 2 

H (p, q) = — muj 4 q 2 



C 



(22) 



where C is an arbitrary constant. The effective distribu- 
tion that will be sampled by the Langevin dynamics can 
be obtained analytically and is 



P e (p, q) oc e 



(23) 



This solution can be normalized only if At < 2/uj. For 
longer time steps, the dynamics is unstable. Although 
the logarithm of the distribution in Eq. ([23]) has an ex- 
pression similar to that of the so-called shadow Hamilto- 
nian for the harmonic oscillator (3ll . |32| , our derivation 
of Eq. (f2"3"| is based on the stationary distribution only 
and does not provide any information on the effective 
trajectory. 

Since H is a state function, it will not exhibit drifts. 
Its square fluctuations can be obtained analytically and 
are equal to 



AH 



1 



1 



2/3 2 



1 



(24) 



For a comparison, the fluctuations of the total energy H 
are 1//3 2 . Interestingly, the size of the fluctuations of H 
depends only on the ratio between the time step and the 
period of the oscillator. It is completely independent of 
the value of the friction. 

The properties of a iV-dimensional oscillator can be 
easily obtained by recalling that, when the dynamics is 
projected on the eigenmodes of the oscillator, the coor- 
dinates evolve independently of each other. Assuming a 
spectrum of N frequencies u>i, the fluctuations of H are 



AH' 



=—y- 



i 



1 



At 4 



N 



(25) 



The last approximation holds when At is much smaller 
than the period of the fastest mode. In this case, it is 
interesting to note that to have rigorously the same ac- 
curacy, the time step has to be chosen proportional to 

JV-V4. 




FIG. 1: (color online) Effective-energy drifts for different 
choices of the friction coefficient, respectively 7=1 (a), 7=5 
(b) and 7=20 (c), and different choices of the time step At, 
as indicated. The effective energy drifts linearly, and its slope 
is strongly dependent on the time step. All the quantities are 
in Lennard- Jones reduced units. 



B. Lennard- Jones crystal 

In the harmonic oscillator the effective energy reduces 
to a state function and does not exhibit drifts. In this 
sense, the harmonic oscillator cannot be considered as a 
prototype of a real molecular system. In this subsection 
we discuss the application of Langevin sampling and of 
effective-energy monitoring in the context of atomistic 
simulations. We use as a realistic test-case a Lennard- 
Jones solid, close to the melting point. We express all 
the quantities in reduced units [1(| • We simulate a cubic 
box with side 19.06 containing 6912 particles arranged 
according to an fee lattice, which corresponds to a den- 
sity p=0.998. We set the temperature to T=0.667. We 
calculate the forces using a distance cutoff of 3. We com- 
pare simulations performed using different values for the 
time step At and the friction 7. All the simulations were 
performed using a modified version of the DL POLY code 

Hi El. 

In Fig. Q] we show a time scries for the effective energy 
H per particle. The effective energy exhibits a regular 
drift due to the finiteness of the integration time step, 
similarly to the total energy in a microcanonical simula- 
tion. The drift is strongly dependent on the time step, 
and is only slightly affected by the choice of the friction. 
In Fig. fj]we show the values obtained for the average po- 
tential energy per particle and the average pressure for 
different choices of 7 and At, obtained from runs of length 
2500 time units. The values are again rather independent 
of the choice of 7. This is remarkable, considering that 
we are changing the friction over three orders of mag- 
nitude and that we are working also in regimes where 
-fAt is not negligible. This indicates that the errors are 
essentially coming from the integration of the Hamilton 
equations and not from the friction itself. In the third 
panel we show the average slope of the effective energy 
per particle, obtained with a linear fitting. The slope is 
again strongly sensitive to the time step and only slightly 
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FIG. 2: Average values of (a) the potential energy per par- 
ticle, (b) the instantaneous pressure and (c) the slope of the 
effective energy per particle, plotted as functions of the fric- 
tion 7. The calculations are performed with different time 
steps: At = 0.0025 (o), At = 0.005 (x), At = 0.01 (A), 
At = 0.015 (□) and At = 0.02 (V). All the quantities are in 
Lennard- Jones reduced units. 
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FIG. 4: Normalized autocorrelation function of (a) the po- 
tential energy and (b) the instantaneous pressure, for different 
choices of the friction coefficient 7, as indicated. The fastest 
decorrelation is observed when 7 is set to an optimal value of 
20. All the quantities are in Lennard- Jones reduced units. 



one wishes to calculate. In particular, to optimize the 
efficiency, the autocorrelation time of the quantity of in- 
terest has to be as short as possible [19|. In Fig. [4] we show 
the autocorrelation function of the total potential energy 
and of the instantaneous pressure, for different choices of 
the friction 7, using a time step At = 0.0025. For both 
the considered quantities, the optimal choice for 7 is 20. 
This rule is far from general, but illustrates clearly the 
fact that too high a friction can spoil the quality of the 
sampling since it hinders particle motion. 



IV. CONCLUSION 



FIG. 3: Average value of (a) the potential energy per par- 
ticle and (b) instantaneous pressure, plotted as functions of 
the slope in the effective energy drift per particle. The calcu- 
lations are performed with different time steps: At — 0.0025 
(o), At = 0.005 (x), At = 0.01 (A), At = 0.015 (□) and 
At = 0.02 (V). All the quantities are in Lennard- Jones re- 
duced units. 



dependent on the friction. 

To stress the fact that the effective-energy slope is a 
correct indicator of the integration errors, we show the 
same data in Fig. [3J There, we plot the value of the 
observable quantity as a function of the slope in the 
effective-energy drift. The two quantities are highly cor- 
related, indicating that the effective-energy slope gives a 
realistic estimate of the errors due to the finite time step. 

Up to now we have discussed the sampling accuracy, 
which measures the systematic errors due to the finite- 
time-step integration. In practical applications also the 
sampling efficiency, which measures the statistical error 
due to the finite length of the simulation, is relevant. The 
sampling efficiency depends on which specific observable 



In conclusion, we have studied the properties of a very 
simple integrator for the Langevin equation, derived em- 
ploying the Trotter scheme commonly use in the deriva- 
tion of multiple-time-step integrators. Moreover, we have 
used the concept of effective energy, introduced in a pre- 
vious paper, to asses on the fly the accuracy of this 
integrator in practical cases, ranging from simple one- 
dimensional oscillators to a Lennard- Jones crystal. Fi- 
nally, we have shown how to monitor the effective en- 
ergy in practice. Our formalism can be easily generalized 
to the description of other stochastic dynamics, such as 
dissipative-particle dynamics [29(. 



APPENDIX A: HIGH FRICTION LIMIT 

When 7 — > 00 the integrator in Eq. (| L 2[) can be rewrit- 
ten in terms of the position only: 



q(t + At) = q(t) + f(q(t))— + AU—R (Al) 
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Now, defining D — this equation becomes 

q(t + At) = q(t) + D/3f(q{t))At + V2DAtR (A2) 

which is exactly the Euler integrator for the overdamped 
Langevin equation 

dq(t) = Df3f(q(t)) dt + V2D dW(t). (A3) 

It is worth noting that the increment of H as defined in 
Eq. (|20|) does not depend on 7, and is still valid. In terms 
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